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Abstract 

We give a self-contained exposition of some mathematical aspects of the Mueller-Stokes formal- 
ism. In the first part we review some basic notions of linear algebra and establish a proper notation. 
In the second part we introduce the Mueller-Stokes formalism and derive some useful mathematical 
relation between physical quantities. Finally a useful decomposition theorem is reviewed. 
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I. INTRODUCTION AND NOTATION 



In these notes we have collected some mathematical results that are not easy to find in the 
hterature. We assume that the reader is already knowledgeable about the Mueller-Stokes 
formalism. All the results presented here can be found in the following references: 

[ 1 ] R. W. Schmieder, "Stokes- Algebra Formalism", J. Opt. Soc. Am. 59, 297-302 (1969). 

[2 ] S. R. Cloude, "Group theory and polarisation algebra", Optik 75, 26-36 (1986). 

[ 3 ] K. Kim, L. Mandel, and E. Wolf, "Relationship between Jones and Mueller matrices for 
random media", J. Opt. Soc. A 4, 433-437 (1987). 

[4 ] S. R. Cloude, "Conditions for the physical realisability of matrix operators in polarimetry" , in 
Polarization Considerations for Optical Systems II, R. A. Chipman ed., Proc. Soc. Photo- 
Opt, lustrum. Eng. 1166, 177-185 (1989). 

[5 ] S. R. Cloude, "Lie Groups in Electromagnetic Wave Propagation and Scattering", Journal of 
Electromagnetic Waves and Applications 6, 947-974 (1992). 

[6 ] D. G. M. Anderson and R. Barakat, "Necessary and sufficient conditions for a Mueller matrix 
to be derivable from a Jones matrix", J. Opt. Soc. A 11, 2305-2319 (1994). 

Different authors use different notations which makes difficult to recognize the same result 
appearing on different papers. For this reason we have tried to simplify and unify the 
notation by adopting one which seems (at least to us) to be the closest to the physics 
(especially the quantum physics) of the problem. For example, we have not adopted the 
awkward "optical" notation for the Pauli matrices 
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but we have adopted the standard "quantum" notation 




Of course, as a consequence of this choice, also the Stokes parameters defined in these 
notes are different from the standard "optical" one. If with E = Xx + Yy we denote the 
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electric field of an homogeneous plane wave propagating along the axis z, then our Stokes 
parameters {Sq, Si, S2, S3} are defined as 



(3) 



where the last three columns display the traditional {{I,Q,U,V}), the "Born- Wolf 
{{S^'^,Sf^,Sf'^,Sf'^}), and the "van de Hulst"^ ({Sq", 5?, ^2", ^l^}) definitions of the 
Stokes parameters, respectively. It is curious to notice that this change of notation was 
already suggested in the sixties [1] but it was unadopted. The last three Stokes parameters 
form a Cartesian coordinate system on the Poincare sphere (Fig. 1). This change of notation 
for the Stokes parameters also causes a change in the definition of the Mueller matrix. For 
example, if we write the Stokes vectors in ours and in the "van de Hulst" notation as 
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then, from Eq. (jS)) it is easy to see that S and S^ are related by a unitary matrix Q, 

S = QS^, 
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where 
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Now, let us consider a linear optical process described in the two different notations as 
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^ M. Bom, E. Wolf, Principles of optics, 7th ed., (Cambridge University Press, Cambridge, 1999). 
^ H. C. van de Hulst, Light Scattering by Small Particles, (Dover Publications, Inc., New York, 1981). 
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Then, it is easy to calculate 

= M5^n, 

from which it follows 

M = QM^Q-\ (9) 

This is the sought relation between our definition of Mueller matrix, and the optical one. 

These notes aim to be, from mathematical point of view, as self-contained as possible; all 
formulae are derived, the only omitted derivations are the ones which reduce to an explicit 
calculation. For example the formula (Ji(J2 = ia^ cannot be "demonstrated", it must be 
checked by explicit calculation from the definition of the Pauli matrices. However, all the 
omitted exphcit calculations can be easily done in few seconds with a computer program 
like Mathematica. 

As we already said, these notes focus on the mathematical aspects of the Mueller formal- 
ism, so no emphasis is given to any physical process. For this reason in the first part of this 
script we almost exclusively deal with the case of deterministic (or Mueller- Jones) Mueller 
matrices which requires the knowledge of the same amount of linear algebra results as the 
more general case. However, all formulae derived here can be straightforwardly extended to 
the case of non-deterministic Mueller matrices. 



A. Notation 

A few words about the notation. We use three different kind of indices: Latin, Greek 
and Calligraphic. Latin indices k, . . . run from to 1 and label the components of 2 x 2 
matrices and 2-D vectors. Greek indices /x, z/, a, . . . run from to 3 and label the components 
of 4 X 4 matrices and 4-D vectors. Finally, Galligraphic indices A,B,C,... run from to 15 
and label the components of 16-D vectors. In these notes the Einstein summation convention 
is used extensively, that is the sum on repeated indices (Latin, Greek and Galligraphic) is 
understood. For example 
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a,^ = ^ij.ubu <^ = ^ A^„b^. (10) 

i/=0 
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We often use the direct product of two matrices A and B, indicated with the symbol 

C = A®B. (11) 
For this kind of matrix product, the standard convention for the indices is the following: 

Cikjl = aijhl- (12) 

It worths to note the order of the indices j and k in both sides of this equation; it will play 
an important role in these notes. 

II. MATRIX BASES 

In this section we study two different ways to represent 2x2 matrices and the relations 
between different representations. 



A. The Standard Basis 

Let A e C^^^ denotes a 2 x 2 complex-valued matrix defined in terms of its elements 

[A]ij = aij, = 0, 1) as 

aoo ctoi 

Oio Oil 

Any 2x2 matrix can be put in one-to-one correspondence with a complex 4- vector a e 
by writing 
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(15) 

This rule is very simple and can be easily extended to n x n matrices by defining 

O-ij = (^ni+j, (16) 
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for i,j = 0, . . . ,n — 1. This rule is so important that in the remaining part of these notes 
we shall refer to as the "Rule". At this point it is important to notice that when we write 
a vector a as in Eq. (fT^. we are implicitly assuming that its components a^, /i = 0, . . . , 3 
are referred to the so-called standard basis in R^, that is 
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Analogously Eq. (fT^ naturally suggests the possibility to write 
A = 
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(18) 



= a^e(^), (/i = 0, 1,2,3), 

where summation on repeated indices is understood and the basis matrices e(^) G M^^^ are 



defined as 
e(o) 



(19) 
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Then the numbers {a^}, that we have found by using the Rule Eq. (fT^ . appear to be the 
components of the matrix A with respect to the basis {e(^)}. In order to demonstrate this, 
it is necessary to define a norm in the vector space C^^^ of the complex 2x2 matrices. It 
is possible to introduce a norm in C^^^ by defining the scalar product {A, B} between two 
matrices A and B as 

{A,B} = Ti{A^B} 

= a*jb,j, = 0, 1), 

where summation on repeated indices is again understood and A^ denotes the Hermitian- 
conjugate of A; that is A^ = (A^)* = (A*)'^, where A* and A'^ are the complex conjugate and 
the transpose of A respectively. Moreover, since {A, B}* = {a*jbij)* = aijb*^ = Tr{i?"''yl}, the 
result {A, B}* = {B, A} follows. By explicit calculation, one can see that the basis vectors 
cire orthonormal with respect to that norm: 



(20) 



e^} = Tr{ef .e(^)} 
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where n,u = 0, . . . , 3 and e|^^ = e^-j follows from Eq. (fT^ . Now, having introduced the 
norm Eq. (j20|) . it is easy to calculate the components of the matrix A with respect to the 
basis {e(;x)} as 

= a^{e(^^),e(u)} (22) 

where /x, = 0, 1, 2, 3. 

Then we have shown that it is possible to associate with any matrix A G C^^^ a vector 
a G and there are two different (but equivalent) ways to calculate a: we can either use 
the Rule given in Eq. (fT^ 



(^11=2 



(^,J = 0,1) 



or calculate explicitly 



ttf, = {e(^),A} 
= Tr{6f^)A}, 



(/i = 0,l,2,3). 



(23) 



(24) 



Until now, A was left arbitrary, therefore Eq. (j24j) holds for any 2x2 matrix. If A 
coincides with one of the basis matrices then Eq. gives the components [e(a)]^ 

(/i = 0, . . . , 3), of the matrix with respect to the basis {e(^j)}: 

= {e(M)>e(a)} 



Tr{eJ,)e{a)} 



■'fia-i 



where a,fi = 0,1,2,3. Therefore we can build the basis 4-vectors e(o) G 
the basis matrices e^a) as 
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It is trivial to calculate from Eq. that 
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that is {e(^)} is simply the standard basis in M^. In summary, we have shown that there is a 
one-to-one correspondence between the standard basis {e(^) G M^^^} and the standard basis 



B. The Pauli Basis 



Another basis commonly used in physics is the so called Pauli basis constituted by the 
2x2 identity matrix and the three Pauli matrices. Here we use a normalized version of the 
Pauli matrices defined as 
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An explicit calculation shows that they satisfy the following multiplication table: 



V2cr(At)0"(i/) 


^(0) 


^(1) 


(T{2) 


^(3) 


^(0) 




^(1) 


^(2) 


^(3) 


^(1) 






ifT(3) 


-ic^(2) 


^(2) 


^(2) 






ifX(i) 


^(3) 


^(3) 


ia(2) 


-ia(i) 


^(0) 



Moreover, again an explicit calculation shows that these matrices are orthonormal with 
respect to the norm defined in Eq. ((20)): 



(29) 



where /i, z/ = 0, 1, 2, 3. The Pauli basis is complete; in order to show this we have to calculate 
the components [c'"(^)]a of the matrices (J(^) with respect to the basis {e(^)} in the way we 
learned in the previous subsection (see Eq. with A = cr(^)): 



a- 



Tr{ef„)^M} 
A. 



(30) 



8 



A 



1 

7! 



(31) 



where /i, a = 0, . . . , 3, and in the last hne we have defined the 4x4 transformation matrix 
A in terms of its elements [Aj^^ = A^^ = Tr{e^^cr(^)}. An explicit calculation shows that 

/l l\ 
1 -i 
1 i 

\i -ly 

where Eqs. fll9l28p have been used. In the previous subsection we shown how to build the 
basis vectors {e(^a)} in associated to the basis matrices {£(«)} in M^^^. Analogously, we 
can now build the basis vectors in associated to basis matrices {cr(^)} in C^^^. To 

this end, for a given /i we define the four components (a = 0, . . . , 3) of the vector 

S(p), as = [cr(^)]a, that is 
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(32) 



where /i = 0, . . . ,3. From Eq. (jHzj) is clear that the yU-th column of the matrix A is made 
of the components of the 4- vector Alternatively we can find the vectors by 

using the Rule: [s(^)]a=2i+j = For example, for /i = 2 we have 



^(2) 
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from which we deduce that 
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(34) 



From Eq. ()3ip it is easy to check by explicit calculation that A is unitary, that is A^A = 
AA^^ = I4, where I4 is the 4x4 identity matrix. In terms of the components with respect to 
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the basis {e(^)} the relation I4 = AA''' becomes 

3 

3 

= E ^-^^h (35) 

fj.=0 

3 

The first and the last line of Eq. (jH^ give us the completeness relation (also called resolution 
of the identity) we are seeking: 

3 

X1K')]"Km)]/3 = "^"Z?- (36) 

By using Eq. ()30p this relation can be written in less involved form as 

3 3 



11=0 fi=0 



(37) 



which is easier to understand. It is useful for later purposes to write the completeness 
relation in terms of the "Latin" matrix elements {i,j = 0, 1). To this end we first 

associate the four indices i, j, fc, Z = 0, 1 to the two indices a, /5 = 0, . . . , 3, by using the Rule: 

a = 2i + j, , ^ 

(38) 

P = 2k + l. 

Then, after noticing that Saf3 = S2i+j,2k+i = ^ik^ji, we rewrite Eq. (jH^ as 

3 

^Wit^)]iA(^ii,)]ki = ^ikSji. (39) 

^J.=o 

From the definition in Eq. (j3(J|) . it is obvious that A is the matrix that performs the change 
from the Pauli basis {cr(^)} to the standard basis {e(^)}: 



(40) 
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where /i, = 0, . . . , 3 and summation on repeated indices is understood. Previously we 
learned that to any matrix corresponds a vector, therefore the matrix A also performs 
the change from the basis to the standard basis {e(^)}. In fact, since by definition 

[s{^l)]a = [<^{f,)]a = Aa^, then 

S{fM) = e(^a)[S(f^)]a = e(a)Aa^, (41) 

where a, /i = 0, . . . , 3. This relation can be written in fully matrix form by noticing that if 
U denotes the unitary transformation between the two basis: = f/e(^), then it follows 
that UafM = (e(Q,), ?7e(^)) = (e(Q),S(^)) = Aq,^, where the parentheses symbol {u,v) indicates 
the ordinary Euclidean scalar product in C" 

n— 1 

{u,v) = J2<Va. (42) 

So we have found that U = A and, therefore, 

S(;.)=Ae(^). (43) 

III. THE MUELLER FORMALISM 



Let us consider a doublets of stochastic variables. 



(44) 

which transform under the action of a deterministic optical device, as 

E ^ E' = TE, (45) 

where the 2x2 complex-valued transformation matrix T is known as the Jones matrix 
representing the optical device. We do not make any hypothesis on the nature of the matrix 
T, it can be arbitrary. The quantities Eq,Ei in Eq. ()44j) are complex random variables 
described by a given ensemble. Starting from Eq, Ei we can build the covariance matrix (or 
polarization matrix) J G C^^^ whose elements are defined as 

J,, = {E,E*) (z,j = 0,l), (46) 
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where (■) denotes the ensemble average. Note that this average has nothing to do with any 
random medium, at this stage we are just considering two components of the electromagnetic 
field as two stochastic variables. By definition, J is Hermitian and nonnegative (or, positive 
semidefinite) , that is (x, Jx) > 0, Vx G C^: 

= x*{E,E*)x, 

= i^*E,E*x,) ^^^^ 

= {ix*E,)ix*E,r) 

= {\x*E,f) 

= (|(x,E)|^)>0, 

where i,j = 0, 1 and summation on repeated indices is understood. Moreover, in the third 
line of Eq. (j47p we have exploited the fact that, by hypothesis, the vector components Xi 
are deterministic variables and, therefore, are not affected by the ensemble average. 
As any other 2x2 matrix, J can be written in the basis as 

J = S^a^,) (/i = 0,...,3), (48) 

where the components = Tr{cr(^) J} of the 4- vector S are known as the Stokes parameters 
of the field. Explicitly 

j^^fSo + Ss S,-iS,\ 
V2 ys, + iS2 So-SsJ 
Form the formula above we see that 

Tr J = v^^o, (50) 

while from the definition Eq. ()46p we have 

TiJ={\Eo\') + {\E,\')^I, (51) 

where with / we denoted the total intensity of the beam. By equating Eq. (j^Uj) with Eq. 
fl51|) we obtain our definiton of Sq: 

5„ = ^. (62) 
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Under the transformation T, the polarization matrix J transform as J — J' where, by 
definition, 



4 = ^rE'P 



Tik{EkEi*)T*i (53) 



TikJkiTlj, 



or, in matrix form, 

f = TJTl (54) 

From Eq. (j54j) is clear that the transformed coherency matrix J' is still Hermitian and non- 
negative. In correspondence to the transformation J ^ J' , the Stokes parameters transform 
as — i> S'^ where, by definition, 

= Tr{a(^)J'} 
= Tr{o-(^)TJTt} 

= Tr{(T(^)T5,a(,)Tt} (55) 
= Ti {a (^^)Ta(^)T'i}Su 
= M^,S,, 

where Eq. ()48p has been used in the third line and we have defined the 4x4 Mueller matrix 

M as 

M^, = Tr{a(^)T(T(,)Tt} ^^^^ 
= {or(^),Ta(^)Tt} , 
where /z, = 0, . . . , 3. It is easy to see that M has real elements: 

M;^ = {or(^),TaHTt}* 

= Tr{Ta(,)Ttor(^)} (57) 

= Tr{a(^)T(T(^)Tt} 

= 

where the cyclic property of the trace: Tr{Ai?} = Tt{BA} has been used. En passant 
we may note that if we write the Jones matrix T in the Pauli basis as T = CaO'{a), where 
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Ca = Tr{cr(Q,)T}, then Eq. can be written as 

M^^ = Tr {o-(^)T(T(^)rt} 

} (58) 

= Tr{CT(,.)}, 

where the cychc property of the trace has been used and we have defined the coherency 
matrix C : C/sa = C/jC* and the 16 matrices {r(^,,)} : [r(^j,)]a/3 = Tr {o-(a)cr(^)0-(/3)(T(i,)}. In 
the remaining part of these notes we shall derive again the result in Eq. ()58|1 in two other 
different ways which are perhaps more complex but also more physically clear. 
Note that from Eq. (jSUj) it follows 

Moo = Tr {a(o)Ta(o)Tt} = ^Tr {TT^} , (59) 

therefore, when T is unitary Tr |TT^} = Tr {I2} = 2, which implies 

Moo = 1. (60) 

This is then the "natural" normalization of M. 



A. Prom the M matrix to the H matrix 



The Mueller matrix M has not, in general, any particular symmetry property. However 
it is possible to extract from it an Hermitian matrix H in the way we are going to show. 
Let us start by writing M in component form as 

Mf,^ = Tr {o-(^)To-(,,)Tl'} 

= TnpT,l^g[a(^^)]rnnW(i,)]pq (61) 

(T ® T \nrn,pq\^ {^)\mn\f^ {u)\pq 
— Fnm,pq\(^ {^)\mn\(^ {i>)\pq 

where we have defined the matrix F G C^^^ as 

F = T®T\ (62) 
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which contains all the information about the scattering process. From Eq. ()62|) it is clear 
that F is not Hermitian, however, we can extract out of it the Hermitian matrix H by doing 
a partial exchange of the tows (Per[ ]) defined in the following way: 



H = Per[F] ^ H, 



np,mq 



nm,pqi 



(63) 



where the indices p and m have been exchanged. This definition clearly requires the matrices 
H and F to be written with four indices, as if they were generated by a direct product of 
two 2x2 matrices; see, e.g., Eqs. ()lllll2j) . However, this is unnecessary; actually after a 
careful examination of Eq. (jUHj) one can easily convince himself (or herself) that the effect 
of the "Per[ ]" operation on an arbitrary 4x4 matrix can be written explicitly in matrix 
form as 

^ao bo ai bi\ 



Per 



fao bo Co do\ 

fli bi C\ di 

a2 62 C2 d2 

\a3 63 C3 dsj 



Co do C\ d\ 
a2 62 as h 

\C2 d2 C3 4/ 



(64) 



This equation can be considered as a definition of the Per[ ] operation alternative to the one 
given in Eq. The advantage of Eq. (jMjl with respect to Eq. is that it does not 

require the 4x4 matrix to be written as the direct product of two 2x2 sub-matrices, but 
it is applicable to arbitrary matrices. 

The matrix H is Hermitian: this can be easily seen by first writing explicitly F in terms 
of the components Tjj of T 

/ Jnn-tnn inn-tni J ni J nn J m -t ni 



00^00 J-OOJ-oi J-Ol-l-00 

Ti rri^ rri rruf rri rri^: rri rri^if 

oo-'io -'oo-'ii -'oi-'io -'01-' 11 

10^00 ^10^01 ^11^00 ^11^01 

\J 10^10 ^10^11 J 11^10 ^11^11/ 



(65) 
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and then by applying the Per[ ] operation to F to obtain H: 

H = Per[F] 



/ rri rji^ rji /-y-i* rji rji^ rji /T-i^t: \ 

'-^00-^00 -^00-^01 -^00-^10 -^00-^111 

J 01^00 ^01^01 ^01^10 J 01 J 11 

10^00 ^10^01 ^10^10 ^10^11 

VllJoO ^11^01 ^11^10 ^11^11/ 



/Too\ 

Toi 

Tio 
VTii/ 



T~'=f= T~'=f= T~'=f^ T~'=f= 

-'oo -'oi -'lo -'ii 



(66) 



= hh\ 

where the diad hh) is written in terms of the 4-vector h defined as 

/roo\ 



h 



Tio 
\TuJ 

which is just the 4-vector representing T in the basis {e(^)}: 



T = hne 



(67) 



Then, by using Eqs. (j66H67p we can write H in component form as 

from which its Hermitian character is evident. Finally, by combining Eq. (jU^ and 

get 

F = /i^e(^) (g) hle(^„) 
= H^ye^^,^ (g) 6(1,), 

which shows that H is just the representation of F in the basis {e(^) ® ^{v)}- 



(68) 



(69) 



we 



(70) 
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Now we can continue the calculation of M^^, by inserting Eq. ()62p in Eq. (jHH) obtaining 

= [<yl^)]nm{T ®T*)nm,pq[C'{u)\pq (71) 

where in the second line we exploited the fact that the Pauli matrices are Hermitian, so 
[o"(^)]mn = [<^(*^)]nm, and in the last line we used the Rule to define a = 2n+m and /? = 2p+q. 
But since = Aq,^, then 

= Kl^{T®T*UKp, (72) 
= [At(T®T*)A]^,, 

or, in matrix form 

M = K\T ®T*)K. (73) 

This formula is particular relevant because it permits us to define the matrix F even when 
the Mueller matrix is nondeterministic or, equivalently, when is not a Mueller- Jones matrix. 
In fact, by rewriting Eq. (fTSj) as 

M = A^FA, (74) 
it is clear that we can invert it and define, in the general case 

F = AMAl (75) 

In the same spirit we can define H in the general case by starting from the last line of the 
Eq. (j?)T|) which can be rewritten with the help of the Eq. (jH!^ as 



-^fiu Hnp^mq\(^{fj,)]mn\(^{u) 



Hnp,mq [^(/i)]mn 

Hnp,mq['^{lJ,) ® '^(^^^]mq,np (76) 



where we used the Rule to define a = 2n + p and (3 = 2m + q. It is simple to invert this 
equation by using the completeness relation Eq. ()39|) that here we rewrite 

3 

^W{f^)]ijWlf,)]ki = SikSji. (77) 
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Then, by multiplying both sides of Eq. (frB|) per [o'(^)]fej[c"(*y)]j7 and summing on /i and i/, we 
obtain 

Hnp,raq^mk^ni^pj^ql ('^^) 

which is the desired result. This equation can be put in matrix form by noticing that 
which permits us to write 

HijM = M^^((T(^) (g) cr*^))ij,ki- (80) 

This formula is not very appealing because it contains both Latin indices which run from 
to 1, and Greek indices which run from to 3. This problem can be solved by using again 
the Rule to define a = 2i + j and P = 2k + 1. Finally, we can rewrite Eq. (IHU)) as 

Hal3 = M^^{a(^i,) (g) 0-^^))a/3, (81) 

or, in matrix form 

0,3 

if = 5^M^,(a(^)®a^,)). (82) 
We can consider this formula as the definition of H for arbitrary M. 



B. The Coherency matrix C 

The relation between H and M is linear but quite involved, as can be seen by writing 
explicitly H in terms of the components M^,^ of M: 

Hoo = |(Moo + Mo3 + M30 + M33) , 
Hoi = \ (Moi + M31 + iMo2 + iM32) , 

(83) 

^02 = \ (Mio + Mi3 - iM2o - iM23) , 
^^03 = \ {Mn + M22 + iMi2 - iM2i) , 
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Hl2 — 




-M22 


-iMi2 


- iM2i) , 


Hi3 — 


|(Mio 


-Mi3 


-iMao 


+ iM23) , 


H20 — 


i(Mio 


+ M13 


+ iM2o 


+ iM23) , 


H21 — 


1(^11 


-M22 


+ iMi2 


+ iM2i) , 


H22 — 


|(Moo 


+ M03 


-M30 


- M33) , 


H23 = 


1 / r 


-M31 


+ iMo2 


- iM32) , 


H30 = 


|(Mn 


+ M22 


-iMi2 


+ iM2i) , 


H31 — 


HMio 


-Mi3 


+ iM2o 


- 1^23) , 


H32 = 




-M31 


-iMo2 


+ iM32) , 


H33 = 


|(Moo 


-Mo3 


-M30 


+ M33) . 



(84) 



(85) 



(86) 



From this formula we see that 

Tr{H} = 2Moo. (87) 

If we choose the "natural" normalization Mqo = 1, it follows Tt{H} = 2. The matrix H is 
not the only Hermitian matrix we can extract from M, actually there arc infinitely many 
Hermitian matrices generated by M which differ from H hy a, unitary transformation. A 
particularly relevant Hermitian matrix is the Coherency matrix C defined as the represen- 
tation of F in the basis (8) cr(^)}. In order to find this representation, let us first write 
the transformation matrix T in both the bases {cr(^)} and {e(^)} as 

T ^ c^a^^) ^ h^e^^), (/x = 0,l,2,3), (88) 

and then let us calculate 



F = T^T* 

where we have defined the coherency matrix elements as 



(89) 



C^, = c^c:, (/i,z/ = 0,l,2,3). (90) 
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By comparing Eq. ()70p with Eq. ()89|) it appears evident that C and H are different 
representations of the same matrix F, with respect to different bases. Therefore they must 
be related by a unitary transformation: we want to find it. To this end, let us recall that 
if A is an arbitrary 2x2 matrix in which can be represented in the two different bases 
{e(^)} and {cr(^)} as 

-4 = a^e(^) = 6^cr(^), (/i = 0, 1, 2, 3), (91) 
then the expansion coefficients and 6^ are related by the change of basis matrix A as 



= {e(^),A} 



(92) 



or, in more compact form, 

a = Ab. (93) 

In our specific case we find, by using Eq. 



(94) 



= Cfj_ye(^a)Aafi ® e(/3)A^j, 
= Aa^C^yAl^e(^a) ® e{/3) 
= [ACA%f3e(^a) ® e(/3). 

The comparison of Eq. (fTUj) with Eq. reveals that 

H = ACAl (95) 

Finally, we can combine the results in Eq. (jH^ and ()95|) to obtain the sought relation 
between H, C and M: 

C = A^HA 

0,3 

= J]M^,At(a(^)®a^,))A 

0,3 

where we have defined the 16 Hermitian matrices \V[^^y)\ as 

r(^,)=At(a(^)®a^,))A. (97) 
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Note that since the Pauh basis is complete in C^^^, the direct products {<J{p) ® form 
a complete basis in C''^'^. Moreover, since A is unitary, the matrices cr(^) (g) cr*^^ and r(^j,) 
are equivalent, therefore the 16 matrices {Ti^^^^} form a complete basis in C^^*^. From the 
matrix rule 

{A ® B)^ = A^0 B\ (98) 

and the definition Eq. (j^Tjl . it follows that the matrices {r(^,^)} are Hermitian. Moreover, 
with the help of the general matrix rules 

Tr{A}Tr{5} = Tr{y4 ® B}, (99) 

and 

{A® B){C ® D) = AC ® BD, (100) 
we can show that the matrices ^{^u) are also orthonormal: 

Tr{r(^,)r(,^)} = Tr{At(a(^)®a^^))AAt(a(,)®a^^))A} 

= Tr{At(a(^) ® ^H)(^(a) ® ^(/3))A} 

= Tr{AAt(a(^)a(,)®a*^)a^^))} 

= Tr{(7(^,)(7(«) ®a*^)(T*^)} 

= Tr{(T(^)(T(„)}Tr{(T^^)(T*^)} 

where AA"'" = J4 and the cyclic property of the trace have been used. 

Now we use Eq. fjlOip to invert Eq. (jUUI) and express M as function of C. By multiplying 
both members of Eq. by T(^ai3) and by tacking the trace, we obtain 

0,3 

Tr{r(„^)C} = 5^M^,Tr{r(„^)r(,,,)} 



(101) 



0,3 



= M^p, 

that is 

M^, = Tr{r(^,)C}. (103) 
The Eqs. (j96p and (jl03|) can be put in a more compact form by using the Rule for n = 4: 

(H -^4/i + z/ = ^G {0,...,15}. (104) 
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Then we can rewrite 



where 



15 



^=0 



and where 



mA = Tr{r(^)C}, 



(105) 



(106) 



m 



(107) 



is the 16- vector associated to the matrix M in the basis {r(_4)}. 

We conclude this section by writing exphcitly the relation between the matrix elements 
of C and M. An explicit calculation shows that if C is written as 

( Qq + a c — id h + ig z — ij ^ 

c + id + b e + if k — il 

h — ig e — if bo — b m + in 

\i + ij k + il m — in — a J 



C 



(108) 



then M has the following form: 



M 



^ao + bQ c + n h + I i + f \ 

c — n a + b e + j k + g 

h — I e — j a — b m + d 

\ i — f k — gm — daQ — boj 



(109) 



Alternative version 



In the literature can be found another method, more geometrical, to find the matrices 
r(^j,) and the result shown in Eq. pOH|l . In this subsection we expose that method. 

Let X, Y two matrices in C^^^ and let us consider their product Z = XY. With x, y and 
z we denote the 4- vectors associated to X, Y and Z respectively, with respect to the Pauli 
basis: 

X = Xf,a^^) ^ Xf, = Tr{(T(^)X}, 

Y = y^a(^^) ^ = Tr{a(^)r}, (110) 
Z = z^a(^) ^ Zf, = Tr{cr(^)Z}, 
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where /i = 0, . . . , 3 and summation on repeated indices is understood. We want to find a 
formula which expresses z as a function of x and y. To this end let us write 



111) 



= Tr{cr(^)Z} 
= Tr{a(,.)XF} 

= a;„y/3Tr{o-(^)a-(c,)0-(/3)} 

where we have defined the four matrices T(^) G C^^"* as 

= Tr{a(^)a(„)a(/3)}. (112) 
Then we can rewrite Eq. f)lll|) in a compact form as 

= (f*,T(^)y). (113) 
Let us notice that because of the cyclic property of the trace 

Tr{(T(^)(T(„)(T(^)} = Tr{(T(^)(T(^)or(„)} = Tr{(T(„)a(/3)(T(^)}, (114) 

we can write 

We shall exploit this property in a moment. Now, let us consider the special case in which 
y = o-(^u) =^ yf3 = 5pv Then, from Eq. ()113|) follows that 



= [T(^u)UXa. (116) 

where Eq. (|1 1 5|) has been used. Another special case is the transposed one, that is when 
X = (T(^) =^ Xa = Sa-y and 

The previous results can be summarized as follows: 

Z = ^ z = T(,)f ^^^^^ 

Z = a^^)Y ^ z = Tj^^y. 
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Now we are equipped to consider the last, most complicated case Z = cr(^)Tcr(,^), where 
T = CaCr(Q,) is a given 2x2 matrix associated to the vector c, where 

/co\ 



Cl 

Then, by putting Y = Ta(^y) ^ y = T(j,)C, it is easy to see that 
where Eqs. (|118p have been used. To summarize, we can write 

= r(^^)c, 

where we have defined the 16 matrices r(^,^) G C''^'' as 

-p — 'Y'^r 



(119) 



(120) 



:i21) 



(122) 



and the symbol "=" stands for "is represented by". 

Now we want to use these equations to calculate the matrix elements M^^, of the Mueller 
matrix, by using Eq. ()6H1 that here we rewrite: 

M^, = Tr{TV(^)ror(,)}. (123) 

Before doing that, notice that if T = CaCia) = c then T'l' = c*^<7(a) = c*; and notice that if 
A = a and B = b, then 

{A,B} = Ti {A^B} 



a*6^Tr {or(^)or(^)} 



(124) 



= (a,"5). 

Finally, from Eqs. p21H124|) it follows straightforwardly that 

M^^ = Tr {TV(^)T(T(^)} 
= (c,r(^^)c) 

= Ci3a[T(^u)]al3 

= Tr{CT(,.)}, 



(125) 
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which coincides with the result found in Eq. ()103|) . To complete the calculation we have 
to demonstrate that the {r(^j,)} matrices found in Eq. ()122|) coincide with the ones found 
in Eq. (|97p. To this end we calculate the matrix elements in both cases and then compare 
them. Let us start from Eq. (jl22j) to write 

3 



7=0 
3 



^ Tr {o-(^) a(^) (T(«) }Tr {cr(^) (T(^) cr(^) } 

7=0 
3 

^ ] ij W{j)]jk[cr{a)]ki[cr(u)] Im [0-(7) ]mnW{f3)] nl 

=0 

X]k(7)]ifeK7)]mn K/,)]ii[cr(a)]fc4^H]i»n[0-(/3)]riZ- 



7=0 



v7=0 / 

From the completeness Eq. (jH^ we know that 

3 3 

^[0-(^)]ife[cr(^)]mn = ^[(^{'))\jk[(yl^)]nn 
7=0 7=0 

SO that Eq. ()122p becomes 



(126) 



(127) 



;i28) 



W (m) ] «i ['^(a) ] fe« ['^ (f) ] Im [(^ ] nl 
= Wif^)]ijW{a)]ki[cr{u)]lk[cr((3)]jl 

= Wia)]ki W M]ijWm]jlWH]lk 

= Tr{(T(a)(T(^)a(/3)0-(^)}. 

The equality 

[^fM)^M]"/5 = Tr{a(«)a(^)cr(/3)a(,.)}, (129) 

can be also easily checked by explicit calculation. Now we repeat the calculation of [r(^j,)]Q,^ 
starting from Eq. (jUTjl : 

= [A1,.K)® [A].^ (130) 

where we have used Eq. (|3U|) in the last line. Now we use the Rule to pass from the dummy 
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4-D Greek indices to the dummy 2-D Latin indices and write 

^ W{a)]k^[cr{^,)]^j[(T(^P)]jl[a(^^)]lk 
= Tr{(7(c,)(7(^)(7(;3)(7(^)}. 

This complete our demonstration. 

Let us conclude this subsection by calculating explicitly the 4 matrices {T(^)}. First of 
all we notice that in the case in which one of the indices is zero, then we have 

[T(o)]a/3 = -^Tr {C7(„)C7(^)} = '^^"'l^' 

[^mU = ^Tr{a(^)CT(;3)} = -^5^^, (132) 

1 1 
[T(^)]ao = -^Tr{(j(^)(j(a)} = 

In the case in which all indices are different from zero, we use the following well known 
property of the Pauli matrices 

= {^ij'^(o) + ^£iji<^{i)) , (133) 

where i, j, / = 1, 2, 3 and £123 = -£132 = £312 = -£321 = £231 = -£213 = 1 is the completely 
antisymmetric Levi-Civita pseudo-tensor (all the unwritten components are zero); to show 
that 

= Tr {cTiCTjCrfc} 

= -^^ij^ {(^(o)(^{k)} + ^eijiTr {a^i)a^k)} (134) 
i 
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Finally, by collecting all these results, we can write explicitly: 



T(o) = — 



1 



T(2) = — 



1 



/l o\ 
10 
10 

\0 1/ 



/o 1 \ 

-i 

10 

\0 i / 



^(1) = ^ 



1 

^7! 



1 



/o 1 o\ 

10 
i 

\o -i oy 

/o l\ 

i 

-i 

\1 0/ 



(135) 



IV. THE DECOMPOSITION THEOREM 



In this section we show that a given Mueller matrix M can be written as a linear combi- 
nation with positive coefficients of at most 4 Mueller- Jones matrices. Here we do not adopt 
the Einstein summation convention, therefore repeated indices must not be summed. All 
sums will be written explicitly as in the right side of Eq. (fTUj) . 

In the Eq. (jUUI) we have defined the Hermitian matrix C and we have shown its relation 
with H and M. Since C is Hermitian it can be diagonalized. Let U(p)i ip- = 0, . . . , 3) the 
four eigenvectors of C associated with the four real eigenvalues 



CU(a) = KU{a), (a = 0,...,3). 



(136) 



where there is not sum on repeated indices. The eigenvectors of an Hermitian matrix can 
always be chosen orthonormal, so we assume 



(m(«),M(/3)) = (a,/? = 0, . . . ,3). 

By tacking the scalar product of both sides of Eq. (|136|) with M(^), we obtain 



(137) 



(■U(/3), C'U(a)) = \aiU(i3),U(a)) = K5f3a, («, /? = 0, . . . , 3) , 



;i38) 
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If we write explicitly the left side of this equation we get 

0,3 



0,3 

fJ,,U 

0,3 



where we have defined the matrix U as 



U : Ui3a = [u{a)]p, {a,P = 0,...,3). 



The matrix U is unitary by definition: 



/i=0 
3 



(139) 



(140) 



1411 



= $^Ka)JMl«(/3)JM 
Ai=0 

= (M(„),M(/3)) 
= 5al3- 

By comparing Eq. p38|l with Eq. p39j] we immediately obtain 

[U^CU]pa = \Jpa, (a,/3 = 0,...,3), 

or, in matrix form 

U^CU = D, 
where D = diag{Ao, Ai, A2, A3} or, explicitly 

/Ao \ 

Ai 

D = 

A2 

\0 As/ 

Since C is positive semidefinite, all its eigenvalues are nonnegative: A^ > 0, (/i = 0, . . . , 3). 
Moreover, since from Eqs. (j87l95|) follow that 



(142) 



(143) 



(144) 



Tr{C} = Tt{A^HA} = Tt{H} = 2, 



(145) 
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then Ao + Ai + A2 + A3 = 2. Now we want to write C in terms of its eigenvalues; to this end 
we have to invert Eq. p43|) obtaining 



C = UDU\ (146) 



or, in components form 



0,3 



^ (147) 

3 

= J]A^[n(^)]„[u^^)]^. 
If we indicate with = W(^)-u|^-j the 4x4 Hermitian diad whose elements are 



we can rewrite Eq. ()147|) in matrix form as 



(148) 



3 



C = Y,K^i,)- (149) 

//=0 



It is easy to see that the matrices are orthogonal: 



0,3 

0,3 



fJ,,U 

0,3 

'I 

3 

!/=0 At=0 



(150) 



where in the second line we have exploited the fact that the {^^(a)} are Hermitian. We are 
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now close to our goal; let us notice that from Eqs. ()74ll89ll^ follow that 

M = A'^FA 

0,3 



(151) 

0,3 



Now we insert Eq. (jl49|) in Eq. (jl51|) to obtain 



3 0,3 

a=0 M,^ (152) 
a=0 

where we have defined the four Mueller- Jones matrices (a = 0, . . . , 3) as 

0,3 

These matrices are real, in fact 

0,3 

- Vfo 1 r (^^^) 



since both and r(;^^) are Hermitian matrices. Actually we have still to demonstrate 
that the {$(q-)} are Mueller- Jones matrices. To do that we need two simple partial results. 
The first comes from Eq. ()13H) which shows that 

[r(;,r.)]00 = Tr{cr(o)0-(^)0-(o)Cr(i.)} 

= Tr{a(^)(T(,)}/2 (155) 
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The second result we need is the orthonormahty of the {$(„)}: 

0,3 0,3 
0,3 0,3 
0,3 



(156) 



(157) 



It is now straightforward to calculate from Eqs. ()153I155|) 

0,8 

[*(a)]00 = ^[^ia)]fMu[^(fMu)]00 
3 

fi=0 
3 

fi=0 

= 1/2, 

while from Eqs. (jl56|) we get Tr{$^^$(„)} = 1. A necessary and sufficient condition for a 
Mueller matrix M to be a Mueller- Jones matrix is Tt{M'^ M} = (2Moo)^. In our case we 
have 

% = 1, (« = 0,...,3), (158) 
(2[<l>(«)Joo)^ 

therefore the are genuine Mueller- Jones matrices. This step complete the demonstra- 

tion of the decomposition theorem. In the next subsection we shall derive this result once 
more by explicit construction of the matrices {$(q,)}. 



A. A step backward: from M to T 

Now that we learned how to decompose M, we want to make a step backward in order to 
see if it is possible to find such a kind of decomposition for the 2x2 matrix J' introduced in 
Eq. ()54|1 . To this end we seek a different form for the matrices We start by rewriting 
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Eq. (UnSl) with the help Eq. ^ of as 

0,3 
0,3 



11, u 

3 



(159) 



\At=0 iy=0 ) 

= At (T(,)®T(;))A, 
where we have defined the four 2x2 Jones matrices {T(^a)} as 

3 

The resuh in Eq. fll59|) shows once again that the are genuine Mueller- Jones matrices. 
At this point we can rewrite Eq. (jl52p as 

3 

M = Y, AaAt (T(„) ® Tf,)) A, (161) 

and compare it with Eq. (fTSj) . Then it appears that in the general case comprising also 
nondeterministic Mueller matrices the single Jones matrix T must be substituted by the set 
of the four Jones matrices {T(^a)} following the recipe given above. In the same way, if we 
assume a priori that in the general case Eq. (j54|l must be substituted by 

3 

f = J2^c.T^a)JTly (162) 

and rewrite Eqs. (|61I7HI73'|) . we obtain again Eq. Then the decomposition of J' we 

were looking for has been found. Note that since Aq, > 0, we can always rewrite Eq. 
as 

3 

J' = ^\aT(^a)JTl^) 
3 



= E (a/^^w) J (^^(l)) (163) 

3 



where we have defined A(^a) = \/\J^{a)- In quantum optics and quantum information Eq. 
f|163p is known as "Kraus decomposition" . 
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At this point there are two things to be noted. The first is about the normahzation of 
the Jones matrices T(a). In fact, it is easy to see 

3 3 

fi=0 u=0 



3 



f/;„f/,.aTr{(T(^)a(^)} 



IJ,,u=0 
3 



(164) 



At=0 

= [u^uu = 1. 

This resuh may seem surprising because if the T(q,) were unitary, then the resuh would have 
been Tr{T^|^^T(Q,)} = 2. However, surprising or not, this result is correct and consistent with 
the normalization we adopted. The second thing is about trace-preserving processes. A 
Kraus decomposition maintains the trace of the coherency matrix J, if and only if 

3 

E^wAc-)=^2, (165) 

where I2 is the 2x2 identity matrix. Let us see whether this is true or not in our case: 

3 3 

a=0 a=0 

3 0,3 

a=0 ix,v 
0,3 3 

0,3 



(166) 



where Eq. fll47|) has been used. From the definition Eq. (jHOJ, we can write = (cuC*^), 
where the brackets indicate the average with respect to an ensemble that represent a generic 
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medium. Then, Eq. fjl66p can be rewritten as 



3 0,3 
0,3 

0,3 



V- . \ (167) 



11,1/ 

3 3 



where Eq. ()88|1 has been used. The it is clear that for a non- depolarizing medium TT^ = I2 
which imphes 

3 

E^l)A")=^2. (168) 

a=0 

In summary, for any Mueller M we can calculate the associate Hermitian matrix C. 
Then, by diagonalizing C we find its eigenvectors {u{a)} whose components constitutes the 
Jones matrices {T(q,)}. Finally we can find the transformation rule for the covariance matrix 
J as in Eq. (fTT^ . 

A small comment is in order. Until now we have used the matrix C instead of H because 
it is expressed in terms of measurable quantities. However, from computational point of 
view the use of the matrix H reveals to be more advantageous. This can be seen in the 
following manner: let us multiply both sides of Eq. ()136|) by A and exploit the fact that A 
is unitary: 

(ACA'^) Au(a) = A„Au(„) 4^ Hv(a) = XaV(a), {a = 0, ... ,3), (169) 
where Eq. (j95|l has been used and we have written the eigenvectors V(a) of H as 

^(a) = Au(o), (a = 0, ...,3), (170) 
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At this point we can jump directly to Eq. p6U|) to write T(a) in terms of {7(a) as 

3 

fi=0 
3 

(171) 

3 

where Eq. (pn|) has been used. It is clear then, that the representation of T(a) in the basis 
{e(^)} is very simple, being 

T,., = f l"<-'l° l''°'l'V (a = 0.....3). (172) 

\[Via)]2 [V(a)hJ 

which is very advantageous from computational point of view. 



V. MUELLER MATRIX IN THE STANDARD BASIS 

In this Section we introduce a new Mueller matrix Ai defined with respect to the standard 
basis. Let J and J' be the covariance matrices that describe the input and output light beams 
entering and leaving a given optical system, respectively. We assume the system to be a 
linear, passive optical element described by the linear map Ai: 

M : J ^ J' = M[J]. (173) 

The above linear relation can be explicitly written in terms of cartesian components as 

J[, = M,,,kiJki, it,j,k,le {0,1}), (174) 

or, by using the Rule 

j; = A<^,J„ (/X, I. e {0,1,2,3}), (175) 

where /i = 2z + j, u = 2k + I, and = {e^a), J} = Tr{e^^ J} are the components of the 
covariance matrix J with respect to the standard basis {e^a)}, (« = 0, ...,3). Equation 
f|175p is analogous to Eq. (jHSj), the difference being that the former is written with respect 
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to the standard basis, while the latter with respect to the Pauli basis. Then, it is clear that 
Aifj,u is just the Mueller matrix written in the standard basis. This statement can be easily 
proved by calculating 

J/i = {^{fi), J} 

= Trie! ,J| 

^ (/^) ^ (176) 

= A^.^Tr{cr(^)J} 
= A ^ 

where Eq. fl4()j) was used in the third line (in fact, we have just rewritten Eq. ()92j) ). Now, 
if we insert Eq. ()17f)|l for both Ji, and J'^ into Eq. (jl75j) we obtain 

KaS'^ = M,AupSp, (177) 

which reads, in vectorial form 

AS^' = MAS ^ S' = A^MAS. (178) 

Since we know that S' = MS, then from Eq. p78|l it straightforwardly follows the desired 
relation between A4 and M : 

M = A^MA. (179) 



Finally, from Eqs. ()63l74p it follows that 

M = F, H = FeT[M]. (180) 

It is possible to write Ai directly in terms of the matrix elements of H. To this end, let us 
indicate with {E(^^^)} the standard basis in M^^"^ defined as 

[E(nu)]al3 = S^aSul3- (181) 

An explicit calculation shows that 

Per[E(^,)] = e(^) ®e(,). (182) 
However, this equality can also be easily proved in the following way: Let us write 
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where Eq. (1^ has been used. Now we can use the Rule to write a = 2i + j and (3 = 2k + l 
and rewrite Eq. ()183|) as 

= [e(^) ® ^u)]ik,ji (184) 
= [Per[e(^)®e(,)]]^.^^^ 
= [Per [e(^) ® 

where Eq. fj63p has been used. By comparing the first and the last row of Eq. ()184|) we 
obtain 

= Per[e(^) ® e(^)]. (185) 
Since for an arbitrary matrix A G C^^^ the following relations hold 

0,3 

Per[Per[A]] = A, ^ = A^^E^^^), (186) 

a,l3 

then Eq. p82|) follows and, moreover, we can write 

M = PeT[H] 

0,3 

= J^if„^Per[E(„^)] 



a,p 

0,3 



(187) 



^H^f3 (e(„) (g) e(/3)). 



which is just the sought relation. 



A. 7W as a positive map 

In this subsection, we assume that the linear map A4 is a completely positive (CP) map. 
In this case we can write the transformation law of J as a Kraus decomposition: 

3 

J' = Y,Aa)JAly (188) 

a=0 

In the standard basis 

3 

= (189) 

13=0 
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where 

^/3a = [A«)]/3 = Tr{6;^)A")}- (190) 
If we substitute Eq. (|189|) into Eq. (jl88p we obtain 



0=0 

cj=0 V/3=0 / \7=0 / 

5^^,.<(6(,)J4^)) (191) 



13,1 

where we have defined the Hermitian, positive semidefinite 4x4 matrix x as: 

X = AAl (192) 

Now, in order to compare Eq. (|191|) with Eq. (jl74j) we have to write the latter in terms of 
cartesian components as 

/3,7 



where 



/3,7 

X]^/37[e(/3)]ifc[e(7)]ji > Jki 

/3,7 J 
= ■M.ij^klJkh 

Mij^kl = y^X/37[^(/3)]»fc[^(7)]j' 

/3,7 

= X]x/37[e(/3) ® e(7)]ii,fc^ 

/3,7 



(193) 



(194) 



or, in matrix form 



/3,7 

This Equation should be compared with Eq. ()187|1 to write the identity 

H = x. (196) 
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Therefore, we conclude that when is a completely positive map, its associated H matrix 
is positive semidefinite. 

At this point it may be instructive to write explicitly the relation between M. and x (o^' 
H) in terms of their elements. Since 



5, 



then 



Pa 

= X2i+k,2j+h 



or, in matrix form 



X = H 



( -^00,00 -^00,01 -^01,00 -^01,01 

-^00,10 A^oo,ii -^01,10 A^oi.ii 

-^10,00 -^10,01 -^11,00 -^11,01 

\A<io,io A^io,ii -^11,10 A^ii.ii 



(197) 



(198) 



(199) 



As expected, we found again the relation H = Per[A^], as it is clear from a visual inspection 
of Eq. 



VI. CLASSICAL MUELLER MATRICES AND QUANTUM ENTANGLED 
STATES OR: QUANTUM MEASUREMENT OF A CLASSICAL MUELLER MA- 
TRIX 



In this section we deal with the problem of determining the 4x4 density matrix rep- 
resenting a two-photon state, when the photon pair is scattered by a "medium" classically 
describable by a Mueller matrix. Here, with the word "medium" we denote any linear op- 
tical device, either deterministic or random, which scatters the photons. We consider two 
possible configurations: In the first one, a single scatterer interacts with only one of the two 
photons. In the second configuration there are two spatially separated media, each of them 
interacting with a single photon belonging to the photon pair. The relevant literature for 
the problem under consideration is listed below: 
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[ 7 ] A. Peres and D. R. Terno, J. Mod. Opt. 50, 1165 (2003). 

[ 8 ] N. H. Lindner, A. Peres, and D. R. Terno, J. Phys. A 36, L449 (2003). 

[ 9 ] A. Peres and D. R. Terno, Rev. Mod. Phys. 76, 93 (2004). 

[ 10 ] A. Aiello and J. P. Woerdman, Phys. Rev. A 70, 023808 (2004). 

[ 11 ] N. H. Lindner and D. R. Terno, J. Mod. Opt. 52, 1177 (2005). 

As it is in the style of these notes, we shall follow a didactic approach, so all the main 
formulas will be explicitly calculated step by step. 

A. Rewriting the decomposition theorem 

Let us begin by rewriting Eq. ()162p as: 

3 

J' = ^p„^(,)J^f„), (200) 

a=0 

where we have defined 

Pa = 77^, ^(a) = V2M^oTia), (201) 

in such a way that 

3 

5];p„ = l, Tr{5j„)5(„)} = 2Moo, (a = 0,1, 2, 3), (202) 

0=0 

where Eq. fll(i4j) has been used. Now, we exploit the isomorphism between the classical 
covariance matrix J and the quantum density matrix p and make the ansatz that a single 
photon initially prepared in the quantum state p, after the interaction with a medium 
classically described by Eq. ()200|) can be described by the density matrix p' defined as: 

3 

P = "^Po^S^o^^pSl^y (203) 

0=0 

B. Single- and two-photon quantum states 

Let us denote with 

{K)} = {|0),|1)}, (z = 0,l), (204) 
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the basis kets representing two orthogonal hnear polarization states of a photon. These 
states are often indicated as horizontal \H) and vertical \V), respectively. Here we follow 
the convention 

\0) = \H), \1) = \V). (205) 
By definition these states form an orthonormal and complete basis: 



m = Sij, (z,je{0,i}), ^|z)(z| = i. 

As usual, we put them in correspondence with the standard basis in {/(j)}: 

io) = /w= r h ii) = /(i) 



(206) 





(207) 

In a similar manner, the dual basis (i = 0, 1) is associated with {/(|)}: 

(0|=4 = (10), (l|-/;t) = ( 1 ). (208) 

The two-photon polarization standard basis can be built by tacking the direct product 
between single photon states, as follows: 

\a = 2i+ j) = \i) ® \j) = e {0, 1}, a e {0, . . . 3}), (209) 

where the Rule has been used to write a — 2i + j. It is straightforward to show that 

= ((^|®01)(|A:)®|/)) 
= (^1^)010 

= 5ik5ji (210) 

In the literature it is often used the so-called Bell basis {!&(«))} defined as 

|6(,)) = S|q;), (q; = 0,...,3), (211) 

where the unitary operator B is represented with respect to the standard basis by the 

unitary matrix B 



V2 



/ 1 1 \ 

10 0-1 

011 
yo 1 -1 / 



(212) 
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In explicit form we have 





= 1^0)) = 


1 

V2 


(|oo) + |ii)), 


\r) 


= 1^1)) = 


1 

V2 


(|oo)-|ii)), 




= |6(2)) = 


1 

V2 


(|oi) + |io)), 


ir) 


= l&(3)) = 


1 

V2 


(|oi)-|io)), 



(213) 



where the first column displays the most common notation for the Bell states. 

Four single-photon operators {e(a), (a = 0, . . . , 3)} may be formed by tacking the direct 
product between a single-photon bra and a single-photon ket as follows: 

e(,) = |^)(j|, {a = 2i+3- 2,jg{0,1}). (214) 

These operators can be straightforwardly put in a one-to-one correspondence with the ele- 
ments of the standard basis {e(a)} in M^^^: 

e(«) = e(a) = fii^ ® ^5), (a = 2i + j- t,j e {0, 1}), (215) 

where 

[(^{a)]kl = [f{i) ® 

= S,,Sji (216) 

= ^2i+j,2k+l 

where /3 = 2/c + Z, in agreement with Eq. (|^. 
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C. Two-photon density matrix and scattering processes 



An arbitrary two-photon state can be described by a density operator p as 

0,3 

P = Y,V^p\a){[3\ 

a,l3 
0,1 



^ T^ijM\'^3){kl\ 

i,j,k,l 
0,1 

(217) 



i,j,k,l 
0,1 



i,j,k,l 
0,3 

where fi = 2i + k, v = 2j + 1 and 

Vik,ji = T^ijM ^ ■D = Per[V]. (218) 

At this point we can work directly with the matrix representation of the operators and deal 
with the density matrix p corresponding to the operator p: 

0,3 

where Eq. fj215p has been used. Before going ahead, we need to derive two intermediate 
results. The first one is a simple calculation: Because of the completeness of the Pauli basis 
we can always write: 

3 

0"(a)0"(/i)Cr(/3) = Kafi/SuCTju), (220) 

i/=0 

where, by definition 

Kaiil3iy = Tr{(T(a)(T(^)Cr(^)0-(i,)} = [r(^;,)]a/}, (221) 

where Eq. (11311) has been used. Moreover, we note that from the definition ()22H) it imme- 
diately follows that 

[^{f,u)]af3 = Tr{cr(„)Cr(^)Cr(^)(T(^)} 

= Tr{a(/3)CT( } (222) 
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The second result we need is also a simple calculation: First, from Eq. ()2Uip we write 

3 

^{a) ='^2^'^!^^^/^^' (223) 

(3=0 

where, by definition, 

Sa/3 = Tr{(T(/3)S'(a)} 
3 

= V2Moo ^[%)]^Tr{a(^)(T(^.)} 

M=o (224) 

= a/2Moo [w(a)]/3 

and Eqs. (I2()lll(i()ll4()ll have been used. Then, by using Equations (1401221122411 we can 

calculate the following quantity that will be used later: 

3 0,3 
0,3 



7,o,/3,/i 
0,3 



(225) 



0,3 / 3 



a,l3,fj,,i>,T \7=0 



Moreover, from Eqs. ()142pi43|14(i|l it follows that 

3 3 



7=0 7=0 °° 

A. 



7=0 
3 



^Ua^^X^Ulp (226) 



7=0 
0,3 

7:? 

Co/3; 
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therefore we can rewrite Eq. 

3 



as 



0,3 



7=0 Q,/3,/i,i/,T 

0,3 

0,3 

= E^l'^Tr{CT(,^)}A.,e(.) 

fJ,,U,T 

0,3 
3 

= ^[AMAt],^e(,) 

r=0 
3 

= A<rr?e(r), 



(227) 



T=0 



where Eqs. ()103ll79l226p have been used. 

At this point we have collected all the results necessary to calculate explicitly the trans- 
formation law of the density matrix: 

3 



7=0 
0,3 



0,3 



t 



0,3 



v7=0 



0,3 

0,3 



(228) 



Per 



0,3 

Or 

= Fei[MT>], 

where Eq. ()185|) has been used. Let us note that, by definition, from Eq. ()219p it trivially 
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follows that 

0,3 



a,/3 

= Per[V']. 

Finally, by equating Eq. (j228p with Eq. (j229|) . we obtain 

V' = MV 

(230) 

= AMA^V, 

which, when detlD} 7^ 0, can be inverted to give: 

M = A^V'iVy^A. (231) 

This result shows that the knowledge given by a single input quantum state (represented in 
this case by V) is sufficient to uniquely determine the classical Mueller matrix representing 
the scatterer. 

Equation (j23(J|) relates the Cartesian coordinates in the standard basis of the input and 
output density matrices p and p' respectively. However, in the classical Mueller-Stokes 
formalism the observables are referred to the Pauli basis rather than to the standard one. 
To illustrate this point let us consider the density matrices p^ and p^ of two independent 
photons 

3 

P'' = E^<^^W' iF = A,B), (232) 
and let build the corresponding two-photon density matrix p"^^ in the usual way: 

0,1 

a, 13 ^ ' 

0,1 

= Y^a^^^ia) ® 0-(/3), 

a,/3 

where we have defined the 16 two-photon Stokes parameters as: 

D^^ ^ S^Sl (234) 
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From now on we suppress the superscript AB and we seek the relation between the two 4x4 
matrices V and D defined by the following relations: 



0,1 

0,1 

",/3 



(235) 



By using Eq. ()40j) it trivially follows 



0,1 
0,1 



0,1 



(236) 



0,1 



So, we found 



D = KDK^ ^ D = A^VA*, (237) 



where we have used the fact that A^^A* = /4. Finally, by multiplying both sides of Eq. 
from left by A"'" and from right by A* we obtain 

D' = MD. (238) 

This relation is the "quantum-equivalent" to the classical one relating input and output 
Stokes vectors. Then, in the Pauli basis the expression for the Mueller matrix becomes very 
simple: 

M = D'D"\ (239) 

We can use alternatively Eq. (j231|) or Eq. (j239p to determine what classical Mueller 
matrix is necessary to achieve a certain quantum state. For example, suppose that we seek 
a scatterer that produces a Maximally Entangled Mixed State (MEMS) when interacting 
with an individual photon belonging to an entangled pair prepared in the "singlet state", 
namely |6(3)) as given in the last row of Eq. ()213|) . The output MEMS is characterized by 
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the density matrix ^ in the standard basis 

f 9il) 



v 



\7/2 



where 



^(7) 



7/2 \ 

1-2^(7) 



(7(7) / 

7/2, 7 > 2/3, 
1/3 7<2/3, 



(240) 



(241) 



while the input singlet state is described by 

/o 

1/2 



V 



(242) 



M 



(243) 



o\ 

-1/2 
-1/2 1/2 
\0 0/ 

If we substitute Eq. fl240|) and Eq. ()242|1 into Eq. ()231|) we obtain straightforwardly 

/ 1 l-2g{-i)\ 
-7 
7 
V 1-2(7(7) 1-4(7(7)/ 
As a last example, we consider the case of an output Werner state represented by 

/ (1 -p)/4 \ 

(l+p)/4 ~p/2 
-p/2 (l+p)/4 

(l-p)/4/ 

and again a singlet input state. In this case it is easy to see that the required Mueller matrix 
can be written as 



V 



\ 



(244) 



M 



/l o\ 
j9 
p 
\0 p/ 



(245) 



W. J. Munro, D. F. V. James, A. G. White, and P. G. Kwiat, Maximizing the entanglement of two mixed 
qubits, Phys. Rev. A 64 R030302 (2001). 
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D. Multi-mode states 



Until now we dealt with two- and four-dimensional Hilbert spaces, since we considered 
only polarization degrees of freedom of photons. However, photons also posses other degrees 
of freedom that, although apparently irrelevant, may play an important role. In this subsec- 
tion we consider photons as physical systems with many degrees of freedom, including the 
polarization ones that will be regarded as the relevant ones. 

Let us consider a finite-dimensional "bare bones" version of the electromagnetic field. It 
consists of 2N independent one-dimensional harmonic oscillators each of them characterized 
by two quantum numbers: the "mode" number n e {0, 1, . . . , — 1} and the "polarization" 
number a e {0, 1}. For a given n the two oscillators labelled by the pairs {n, a — 0} and 
{n, a = 1} "oscillate" along two mutually orthogonal directions fixed by the two (possibly 
complex) unit vectors e„o and e^i, respectively: 

(e„a, enp) = 5^p, (a, (5 e {0, 1}). (246) 

A third unit vector e„3 orthogonal to the other two remains automatically fixed by the 
relation 

^n2 — CnO X ^nl- (247) 

It is important to note that in the theory there is not a third harmonic oscillator labelled 
by {n,a = 2} that oscillates along e„2- However, from a geometrical point of view the 
introduction of e„2 is necessary to write the resolution of the identity in a 3-dimensional 
space as 

2 

E ^"i^ii = ^3, (248) 

i=0 

where I3 is the 3x3 identity matrix. A set of N projection matrices {Vn} (and the com- 
plementary ones {Qn}) projecting onto the physical directions of oscillation of the system, 
can be easily build as 

1 

= S'"°'^' (249) 

and Pn + Qn = -^3 • Each harmonic oscillator is characterized by its annihilation and creation 

operators a„a and aj^^^ respectively, that satisfy the canonical commutation rules: 



= ^nrJ>ocP- (250) 
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The Hamiltonian of the system is just the sum of the Hamiltonians of the 2N harmonic 
oscillators: 

^ Af-l 1 

^ = 2 XlXl ^" (aL^na + a-naaia) , (251) 

n=0 a=Q 

where h = 1 and Un > 0. The single-particle states are built from the vacuum state 

|0) in the usual way 

\na) = aL|0). (252) 
Finally, the resolution of the identity can be written as 

I = lo + Ii + • • • 

^^'^ >^ (253) 

= |0)(0| + / ^ / ^ \na){na\ + y ^{multiparticle states}. 

n=0 a=0 

Now that our system is well defined, we try to build a Positive Operator Valued Measure 
(POVM) in order to determine the relevant density matrix pertaining to the relevant po- 
larization degrees of freedom. Let {/i} denotes an orthonormal and complete basis in C^: 

2 

(/i>/j) = '^ij, Y.f^fi=^-^^ (i,jG {0,1,2}). (254) 

i=0 

By using Eq. ()248|) for each mode n we can write 

fi = h- fi 

2 

= ^ ^nj^lij ■ f i 

(255) 



j=o 

2 



j=0 
2 

j=0 

where = (e„j, fi). Then, we define the physical vectors f„i associated to the mode n as 

f ni T^n fi 

1 

= ^ ] ^na^lia ' f i 

(256) 



a=0 
1 



^ ^ ^na (^na; ^i) 

0=0 
1 

= Ve F 



na\i 

a=0 
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(257) 



These vectors are not of unit length nor mutually orthogonal: 

= (/i,^n/j) 

where, by definition, Vnn > 0. 

Now we are ready to write the single-mode operator F„i acting on the physical states of 
the system as 

^ I ^ ^ / , „ „ =F i^nij ^na) ^nai (fni) ^ni) 7^ 0; 

^ ni— \ q;=o V V "i' nij 

(f„i,f„i) = 0. 



(258) 



Then, we can use this operator to build the multi-mode Hermitian positive semidefinite 
"intensity" operator as 



JV-l 



n=0 

AT-l 0,1 

EE 

ri=0 

JV-l 0,1 



^^(fni) fni) 



t 



■\/(fni) fni) 



(fni, Cn/s) ^n/? 



^ ^ ^ ^ (fni, ^na) (fni, ^n/j) ^na^' 



n=0 a,/3 
AT-l 0,1 



^ ^ ^ ^ (Cna, fni) (fni, ^n/?) 



n/3 



in/3 



(259) 



n=0 q:,/3 
JV-l 0,1 

~ ^ ^ ^ ^ (Cna, /ni) (^ni, ^n/j) t*na^"i8' 

n=0 a,/3 

where the last step trivially follows from the fact that Vn^na — ^na and, therefore, 

(^na, fni) {T^n^nai fni) 

(^na, "Pnfni) 
— ( ^nrv. 1 fni ) • 



(260) 



At this point it is easy to see that the three operators {Fq, ^,-^2} form a POVM in the 
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one-particle space: 



2 



i=0 

iV-1 0,1 2 




(261) 



n=0 i=0 

N-1 0,1 




n=0 a=0 



= N, 



where N is the particle-number operator and Eqs. ()246|) and ()254p have been used. 
E. Reconstruction of the density matrix 

Let TZ = {x,y,z} be an orthonormal Cartesian coordinate system in and let U, V 
and W three mutually unbiased bases for defined as 

U = {uo,Ui,U2} = {x,y,z}, 



From a physical point of view, these three bases correspond to the three pairs of mutually 
orthogonal polarization directions {U, V, W: linear horizontal-vertical, linear 45°-135°, and 
circular right-left, respectively), selected by a polarizer whose planar surface is orthogonal 
to z. We want to calculate the Stokes parameter of a beam of light (either classical or 
quantum). To this end, let us imagine to repeat the construction of the POVM outlined 
in the previous section for each of the basis set U, V and W, thus obtaining three different 
POVMs denoted with Ui, Vi and Wi, respectively. For example, if in Eq. ()259p we substitute 
fni with Uni, we obtain 




(262) 



N-1 0,1 





(263) 



n=0 a,l3 
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In exactly the same manner we may obtain Vi and Wi. As a subsequent step we introduce, 
in analogy with classical optics, four Hermitian "Stokes" operators defined as follows: 

1 



(264) 



^(2) = i= (#0 - m 

For sake of clarity, we introduce the six operators {Ex}, {X = 0, . . . ,5) defined as 







' Uo' 


El 




Ui 


E2 




Vo 


Ez 




Vi 






Wo 


A. 




Wi 



in such a way that we can rewrite Eq. (j264j) in a compact form as 



S, 



(A) 



Y^PaxEx, (^g{0,...,3}). 



x=o 



where we have defined the 4x6 matrix P as 

/ 1 1 \ 



P 



1 



V2 



1-10 
1-1 
\ 1 -1 / 

and PP^ = I4. It is instructive to write explicitly the operators {S'(^)}: 



(265) 



(266) 



(267) 



(A) 



PaxEx 



5 



N-1 0,1 

Pax ^ ^ (^na, $,nx) {$,nx, ^n^) (^na^^nlS, 
X=0 n=0 o,/3 



(268) 



where ^ = ${X) G {u, v, w} and x = x{X) G {0, 1}. Then, we can rewrite Eq. ()268|1 as 



N-l 0,1 
n=0 a,P 



X=0 



(269) 
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where an explicit calculation shows that 



t 

nx 



X=0 



\ 



h-A)]oo [^(-^)]oi 

[^(-4)] 10 




(270) 



(A), 



where {<7{a)}: {A E {0, . . . , 3}) are the 2x2 Pauli matrices, and Eqs. ()262|) and ()267p have 
been used. Finally we can write 

N-l 0,1 

n=0 a,l3 ,cy^^ \ 

N-l 0,1 l^'-LJ 

= y^^na, Cr(A)^nf3)aiaanf3, 

n=0 a,f3 

where with G we have denoted the restriction of to a two-dimensional subspace: 



(272) 



Of course, the two-dimensional vectors {Sna} are not unit length nor mutually orthogonal. 
Now we can use this result to calculate 

N-l 0,1 



[niy 



S{A)\mn) = ^^(£pc.,(^(^)£p/3)(ni/|aj„ap^|m/i) 

p=0 a,/3 
AT-l 0,1 



^ ^ (£pa , o-(^) £p/3) (0 1 a„i.aj„ap/3am/, 1 0) 

p=0 a,/3 

^nm(^ni^) '^(yl)^n/j)' 



(273) 



At this point we have all the ingredients necessary to calculate the expectation value (5'(_4)) 
with respect to the generic state described by p: 

0,N-1 0,1 



m,n fi,u 



nv\. 



(274) 
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Then 



0,N-1 0,1 



m,,n fj,,u 
0,N-1 0,1 



Y Y PmfM,nu{nu\S(^A)\miJ,) 

m,n fi,u 
N-1 0,1 

n=0 fj,,u 
N-1 



n=0 



(275) 



(276) 



where we have defined the 2x2 single-mode matrices Vn and <Jn(A) as: 

[(^niA)]aP = {^na,Cr{A)^n(3), 

and a, P E {0, 1}. 

In a paraxial regime of propagation there is a "dominant" mode of the field, say n = no, 
and one can assume 



(£„a, (T(^)£„/3) = {Snoa,0'(_A)£nop), Wn E {0, . . . , N — 1} . 



(277) 



Since we always have the freedom to choose our reference frame, in this case it is convenient 
to choose the two polarization vectors {e„gQ,} associated to the mode uq in such a way that: 



From the definition ()283|) it trivially follows 



e, 



noi = y 



VV 



(278) 



^nnO 



^nnl 



(279) 



which implies 



(,^noa,(^iA)£nol3) - [<^{A)]^f 



(280) 
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Then, from Eqs. (I275I277I28()|) it follows that 

N-l 

i^i^)) = X^Tr{r'„(T„,(^)} 

n=0 

^^'^ T (281) 

n=0 

where we have defined the 2x2 matrix V as 

N-l N-l 
V=^Vn, or, [V]^p = Pna,nP, (282) 

n=0 ra=0 

which coincides with the naive definition of reduced density matrix. 
F. The relevant density matrix 

At this point we are ready to calculate the relevant density operator for the polarization 
degrees of freedom. However, before doing so, let us apply the formulas written above to 
the simple case in which a single mode of the field, say again tiq, is excited. In this case, by 
definition 

Pin^,nu ^mno^nnoPno^,no,u — ^ mno^ nno P ^ (283) 

and if we substitute Eq. (|283|) into Eq. (j275|) we obtain 

0,1 

P = Yplu\n^P)M- (284) 
From Eq. (j275p we can easily calculate 



(^(^)) = Tr{pS(^)} 



0,1 



0,1 

Yplu{noJ^\S{A)Wp) 



(285) 



11,1/ 

0,1 



By substituting Eq. ()280|) into Eq. ()285p we immediately obtain 

0,1 

fJ,,U 

= Tr{pV)}. 
56 



(286) 



since, by definition, Tr{p°} = Tr{/3} = 1, then (5'(o)) = 1/v^ and, after a simple calculation, 
we obtain explicitly 

^^pO^^f (%> + ('^(3)> ('^(l)>-^('5(2)> I 

for later purposes it is useful to define the four scaled real parameters {s^} as 

s^ = i=i§4' (^e{0,...,3}), (288) 
V2 (6(0)) 



p° = ^ + ^^-"^ 1. (289) 



and rewrite Eq. p87|l as 



V2 \ s^ 



+ IS2 So - S3 



Now we go back to the general multi-mode case. Let us suppose that we have measured 
or calculated the four values {(S'(^))}. Two cases are possible: Either 

3 3 
B=l B=l 

If the first case occurs, then it is not possible to calculate a relevant density matrix for 
the system (this may happen because of unwanted experimental errors). Vice versa, if the 
second case occur, then can be obtained straightforwardly. This result can be achieved 
in three steps: First, by using Eq. ()288|1 we calculated the four scaled parameters {s^}; 
second, we introduce the four relevant observables {(3"^^)} such that 

Tr{pX^)}=s^, (291) 

where, inspired by Eq. ()27H) . we make the ansatz: 

^iA) = ^{A), (^G{0,...3}). (292) 

Finally, we calculate the less biased by using the maximum entropy criterion which leads 

to 4 



p^ = p^ = exp 



6=1 



(293) 



R. Balian, Incomplete descriptions and relevant entropies, Am. J. Phys. 67 (12), 1078 (1999). 
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where the normahzation term \1/ ensure the condition Trp^ = 1, and each of the three 
lagrange multiphers 7^, {B = 1,2,3) is associated with each constraint Eq. ()29H) . After a 
straightforward calculation it follows that 

3 

= V^ = J2sAcriA)- (294) 

This is our resuh, now we seek a relation between input and output relevant density operator 
in a multi-mode scattering process. 

G. Input and output relations 

Let us consider now a generic linear scattering process that transform the input single- 
photon density operator into the output single-photon density operator p°"*: 

r' = 

i 

where 

J2AAl = i. (296) 

i 

The relevant quantities to calculate are the transition amplitudes 

{mfi\Ai\nu) = Ai^f,^{m,n), (297) 

where with Ai{m, n) we have denoted the 2x2 matrix whose elements are Ai^^^{m, n). Then, 
we can rewrite Eq. (j295j) as 

o,Af-i 0,1 

pZbp = E E E ^)pi,n.4./3(^' b). (298) 

m,n fijU i 

From the algebra of the Pauli matrices it is easy to see that for any given pair of modes 
{m, n} one can always write 

3 

= ^5^(m,n)[cr(^)]^j„ (299) 

where we have defined 

0,1 

Pmfj,,nu Mut., iAe{0,...,3}). (300) 

tJ.,U 
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If we use Eq. (fM^ in both sides of Eq. we obtain 

0,1 

a,l3 

0,1 0,Af-l 0,1 

a,/3 m,n fi,i/ i 

0,1 0,iV-l 0,1 3 

a,0 m.n fi,u i B=0 ('^^\^ \ 

3 0,Af-l 0,1 I-JUJ-J 

B=Q m,n i a,l3,fi,u 

3 0,Ar-l 

B=0 m,n i 
3 0,Ar-l 

= ^ ^ MM3ia,b,m,n)S'^{m,n), 

B=0 m,n 

where, by analogy with the definition of a classical Mueller matrix, we have defined 

)a(s)4(n,6)}. (302) 

i 

From Eqs. (ETHll and it follows that 

N-l 0,1 



sm)> = EE 

n=0 fi,u 
3 N-l 0,1 

B=0 n=0 
3 N-l 

= X^X]'5B(n,n)Tr{a(B)a„(^)} 

B=0 n=0 
3 Af-1 

= ^ ^ A^(n)5B(n,n), 



(303) 



B=0 n=0 



where 



^ABin) = Tr {cr„,(^)(T(B) } . (304) 
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Note that in the paraxial approximation (see Eq. (j2ZZ|)) ^AB{n) = (J^b- By using the results 
in Eqs. ()301|) and ()303|) we can write 



3 N-1 

B=0 71=0 
3 N-1 3 0,Af-l 

= J^J^^^bHJ^ Mtsc{n,n,p,q)S'^{p,q) (305) 

B=0 n=0 C=0 p,i} 

0,3 0,N-1 

{n, n,p, q)S'(?{p, q). 

B,C p,q,n 

In a realistic experimental configuration, the input beam of light has a well defined polariza- 
tion irrespective of the spatial and temporal coherency properties of the beam itself. This 
means that it is possible to write 

Pmfj,,nu Pmn,iiu — Rmn^fMU^ (306) 

namely 

R^r, (307) 

where R and r are a N x N and a 2 x 2 matrices, respectively, and Tt{R} = 1 = Tr{r}. 
Note that this factorization has been made upon the matrix p representing the operator p 
and not on the operator itself, where it would have been meaningless. With this assumption 
we can write 

0,1 

SciP,q) = YPpl,qpW{C)]/3a 

a,f3 

~ Rii^V^^irir^ 1 ^^^^^ 

a,/3 

pinoin 

— ^pq'^C 1 

where = Tr{r™(7(c)}, and we use this result in Eq. (|305|) to obtain 

0,3 0,N~1 

(n, n,p, q)S'^{p, q) 

B,C n,p,q 
0,3 0,N-1 

= E E ^AB{n)MBc{n,n,p,q)R-;:^^S^^ 

B,C n,p,q 

rO,Ar-l 3 



B y'^^^) 



(n, n,p, q)R^ 



I m 

"pq 

n,p,q B=0 



c=o 

3 

c=o 



(309) 
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where we have defined the effective 4x4 Mueller matrix M as 

0,Af-l 3 

n, n, p, q)R^. 

p,q,n B=0 

From Eqs. ()303|306l308p . it follows that for a paraxial input beam 

3 N-1 



(310) 



S=0 n=0 

3 N-1 

B=0 n=0 

N~l 



(311) 



1 in 

"nn 



Finally, by comparing Eq. 

(^(^)) and {S(^A)) ■ 



n=0 

= 5'^Tr{i?''^} 

cin 

— 

with Eq. ()311|) we found the sought relation between 



(312) 



B=0 



where 



0,Af-l 3 



= E E^^c(ri)McH(n,n,p,g)/?;° 

p,q,n C=0 
0,N-1 3 

p,q,n C=0 i 
0,Af-l 

E Y^"^ {(^{A)Mn,p)(y(B)A\{<i^T^)} 



in 



(313) 



1 in 



p,q,n i 

where the last, approximate equality is valid only in the limit of paraxial detection. 



H. Two-photon scattering 

Let us consider now the case of two photons, say A and B, that are scattered by two 
independent, spatially separated media. We denote with \aa) and \h(3) the single-photon 
basis states for photons A and B respectively, where a, 6 G {0, . . . , — 1} and a, /3 G {0, 1}. 
A two-photon basis state will be indifferently written as 



\aa) ® \h(3) = \aa)\bp) = \aa,h(3) = \AB), 



(314) 
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where A and B are cumulative indices for the pairs of indices (aa) and (6/3), respectively. 
Let denotes the density operator describing the input two-photon state: 

0,N-1 0,1 

p"" = ^ ^ Paa,bi3;a'a',b'f}'\aa,bP){aa',b'P'\ 



a,b OL.j3 
a',b' a' .13' 



(315) 



5^5Zpab,A'B'|^5)(A'5' 



A,B A',B' 

where 

Paa,fe/3;a'a',6'/3' = (a", fc/^lp'^'la'a', = pai3,a' {ab, a'V). (316) 

A linear scattering process due to two independent, spatially separated media, transforms 
the input two-photon density operator p™ into the output two-photon density operator p°"*^: 

where 

Y,AA = i = Y,B,B]. (318) 

The relevant quantities to calculate are the transition amplitudes 

{aa\Ai\a'a') = A,aa'(a, a') = Ai^AA', , . 

mBjlb'P') = B^^pp,(b,b') = Bj,BB', 

where with Ai{a,a') and Bj{b,b') we have denoted the 2x2 matrices whose elements are 
Ai^aa'{o,,o,') and Bj^j3/^/{b,b'), respectively. Then, we can rewrite Eq. ()317|) as 

out 

Paa,bP;a'a',b'f3' 

0,N-1 0,1 



A,aa"{a, a")Bj^i3pM{b, &")P™'a",6"/3";a"'a"',b"'/3"' AU^o'I*^'"' '^')-S],/3"'/3' (^'") 

i,j a",i)" a",/3" 
a"',b"' a'", 13"' 

^«,AA"-Bi,BB"P!4"B",A"'B"'^I,A"'A'-^l,B"'B'- 



A" ,B" 
A"'B"' 



(320) 

From the algebra of the Pauli matrices it is easy to see that if we define the 4x4 matrices 
^(AB) as 

^(AB) = (^{A) ® cr(B), (321) 
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they form a complete (by definition) and orthonormal set of basis matrices in C^^^: 

Tr {S(^B)S(^/f5')} = Tr { ((T(^) ® (T(b)) (cr(^/) (g) cr(B/)) } 

= Tr{(T(^)(T(^/)}Tr {a(B)(T(B/)} (322) 

Then, it is clear that it is always possible to write 

0,3 

(323) 

where we have defined 

0,1 

Syijs{ab, a'b') = Paa,bP;a'a',b'P' ^{AB)\a' 13' ,aP 



a,f3 

0,1 



a'l3',al3 



(324) 



a, 13 
a',f3' 



= Tr {p(a6,a'6')E(^)}. 

If we use Eq. ()323|) in both sides of Eq. ()320p . we obtain, after a lenghty but straightforward 
calculation, 

0,3 0,Af-l 

S°J^{ab,a'b')=J2 Yl M(^b,a'b';a%'\a'%''')]^,^,i,,S%s'ia'b'',a'''b'''), (325) 

A',B' a".b" 
a"',b"' 

where the 16 x 16 matrix M(a6, a'b'; a"b", a"'b"') is defined as: 

M(a6, a'b'- a"b", a"'b"') = M^^^ (aa', a" a'") ® M^^) (66', b"b"'), (326) 

and where, as in Eq. ()302j) . we have defined the 4x4 matrices M^'^^aa' , a" a'") and 
M^^\bb' ,b"b"') as 



M^^l,{ao!,a"a!") = ^ Tr A(a, a")a(^')4(«'", 

Mg]{bb',b"b"') = ^Tr{a(e)5,(6,6")a(so^J(6'",6')|- ^^^^^ 



Now we want to relate the quantities displayed in Eq. ()325|) with quantities that are actually 
measured which, therefore, corresponds to mean values of Hermitian operators. To this end. 
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we calculate step by step the mean value of the two-photon Stokes operator S'(^) ® S(^t3y. 

= E PAB,A'B'\AB){A'B'\ (V) ® S^B)) } 

A,B A',B' 

= E E PAB,A'B'Tr{\AB){A'B'\ ® 5(b)) } 



A,B A',B' 
A.B A',B' 



A,_B A',B' 
0,N-1 0,1 



E E Paa,b/3;a'a',b'l3'{aa'\S(A)\aa){b'P'\S(l3)\bP) 

a,b ct,f3 
a',b' ci',l3' 
,N-1 0,1 

E^ E/ Paa,b/3;a'a',b'l3'Sa'a {^a'a', 0'(A)£aa) h'b {^b'13', 0'(B)£b^) 



(328) 



a,b ct,f3 
a',b' a', ft' 
0,N-1 0,1 



0,iV-l 0,1 

E ^ Paa,bP;aa',bl3' [c^a(^)]„,^ [^b(S)]^//3 

a',/3' 

0,Ar-l 0,1 

a,b a./3 
0,iV-l 

^ Tr{p(a6,a6) (cr„,(^) ®ct6(b))}. 

a,6 



From Eq. ()3U4|) it is easy to see that 



:'="'=o 



(329) 



^rrU) = Tr{ (T/(^)(T(^/)} 
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where / G {a,b}, and JF G {A,B}. By using this result and Eq. 

as 

0,N-1 



we can rewrite Eq. 



[S^A) ® S(^B)) = ^ Tr{p(a6, ab) {aa(A) ® (^b(B)) } 

a,b 

0,N-1 0,3 

= ^ ^ Aj,_^,{a)ABB'ib)Ti{p{ab,ab) (o-(^/) ® o-(b'))} 

0,N~1 0,3 

= ^ A^^/(a)ABB/(6)Tr{p(a6,a6)S(^/B')} 

0,Af-l 0,3 

= A^^/(a)ABB'(6)5^'B'(a6,a6) 

0,N-1 0,3 

AB,A'B' '^A'B'{o,b, ab) 

a,b A',B' 
0,N-1 

= ^ Sj,B{ab,ab), 



(330) 



where the last, approximate equality holds in the paraxial limit where A(a) = /4 = A(6). 
By using Eqs. ()325|) and ()330p it is easy to see that 

0,Af-l 0,3 



{S{A) ® Si^B)) 



out 



E [^(«) ® ^^^)Wa'B' S%'i^h, ab) 



a,b A',B' 
0,Af-l 0,3 



ABA'B' 



a.b A'.B' 
a'.b' A",B" 
a'l.b" 



0,Af-l 0,3 



£1,6 >i'.e' 

a', 6' >t",B" 
a",b" 



0,Af-l 0,3 



a', 6' A" ,B" 
a".b" 



aa, a' a") 



M^^\bb, b'b")]^,^,^^„^„ S%^.{a'b\ a"b"] 



'N-l 



N-1 



Y A(a)M(^) (aa, a! a") ® ^ A(6)M(^) (66, 6'6") 



a=0 



6=0 



J ABA"B" 



xS'Xnis"{a'b',a"b" 

0,N-1 0,3 

= J]] ^ M_AB,A"B"{ab' , a"b")S%,i^„{a'b' , a"b"), 



a',b' A",B" 
a".b" 



(331) 
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where we have defined 



MABA"ts"{ah\a"h") 



N-l 



A(a)M(^) (aa, a' a") ® A(6)M(^) {hh, h'h") 



a=0 



b=0 



AB, 



or, in compact matrix form: 



N-l 



M{a'b',a"b") = ^ A(a)M(^)( 



aa, a' a") 



^A(6)M(^)(66, b'b"). 



a=0 



6=0 



Then, we can rewrite Eq. (|33ip as 

0,Af-l 0,3 



,A'B'{ab\ a"6")5^,g,(a'6', a"6"). 

a'.y A',B' 

a",b" 

When the input state is not hyp erent angled, one can write 



Paa,bf3;a'a',b'f3' — ^ab,a'b' ' al3,a' (3' ^ 



where i?™ and are a A^^ x A^^ and a 4 x 4 matrices, respectively, and Tr{i?™} 
In this straightforward calculation shows that 

0,3 



{S(A)®S(B)) = ^AB,A'B'S'X'B'y 



A'B' 



and 



0,1 



where we have defined 



and 



S'2'l3'{a'b',a"b") - P'^a,b'0-,a"a',b"l3' [^iA'B')]a'(3',afi 



0,1 

a' ,13' 

= -R™b',a"fe"Tr{r"^S(^/e')} 

pin cin 

= Tr{r-S(^,s,)}. 

0,Af-l 
a,b 

— ^AA'^BB', 
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where the last, approximate equahty holds in the paraxial limit only. In this limit, we have 

{S^A)®S^B)y'' = S%, (340) 
and, by combining this result with Eq. ()337p . we obtain 

S%^,{a'b\ a"b") = i?;?,, V') ® S^B')f. (341) 
Finally, we substitute Eq. ()341|) into Eq. ()334|) to obtain 

0,Af-l 0,3 



A',B' 

a",b" 

0,3 

= E ^AB,A'B'{S{A') ® ^(B')) 



a'M A',B' 

a^'M' (342) 



A',B' 

where we have defined the field-dependent two-photon 16 x 16 Mueller matrix M as: 

0,Af-l 



Mas,A'B'= Mj,sA'B'{ah\a!'h")B}^y^,„y, (343) 



a>,b' 
a",b" 



Equation (j342j) is our final result: It represent the linear relation between input and output 
measured quantities. This equation is the two-photon quantum analogue of the classical 
Mueller-Stokes relation. This similarity can be made manifest if we define the 16 two-photon 
Stokes parameters as 

{Si^A)®Si^B)) = S^, (344) 

where we introduced the cumulative index <l> = {AB) G {0, . . . , 15}. Then, Eq. (I342|) can 
be rewritten as 

15 

5T = J]M$$,5i°„ (345) 



*=o 

which is formally equivalent to the classical one. 
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FIG. 1: The Poincare sphere. To each point on the sphere it is possible to associate a definite pure 
polarization state of the light. Moreover, internal points are associate with mixed (or partially 
polarized) states. 
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